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The study of extragalactic sources of high energy radiation via the direct measurement of the 
proton and neutrino fluxes that they are likely to emit is one of the main goals for the future 
observations of the recently developed air showers detectors and neutrino telescopes. In this work 
we discuss the relation between the inclusive proton and neutrino signals from the ensemble of all 
sources in the universe, and the "resolved" signals from the closest and brightest objects. We also 
compare the sensitivities of proton and neutrino telescopes and comment on the relation between 
these two new astronomies. 
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^ . I. INTRODUCTION 

< 

t^J- ■ There is a general consensus that the highest energy cosmic rays (CR) are of extragalactic origin because they 
reach the Earth with an approximately isotropic angular distribution. The magnetic fields of the Milky Way are not 
sufficiently strong and extended to randomize the directions of particles produced by our own Galaxy, and the this 
isotropy of these CR is likely to reflect the large scale homogeneity of the universe. It is natural to expect that the 
extragalactic Ultra High Energy Cosmic Rays (UHECR) are produced in "point-like" astrophysical sources. The 
identification of these sources, and the clarification of the mechanisms that accelerate particles to these very large 
energies is clearly one of the crucial goals of the new large acceptance air shower detectors like the Pierre Auger 
Observatory. 

The trajectories of charged particles is bent by the presence of astrophysical magnetic fields, however at sufficiently 
large rigidity E/Z (with Z the particle electric charge) the magnetic deviations should become sufficiently small to 
allow the direct imaging of sources with cosmic rays. Unfortunately the intensity and structure of the intergalactic 
magnetic field are very poorly known, and therefore the region (in source distance and particle energy) where CR 
' astronomy is possible remains very uncertain. It is possible (and there are observational hints) that the CR above 
£T) . E > few x 10 19 eV are already propagating in quasi-linear mode for source distances of order 100 Mpc or more, 
however the identification of the sources remains difficult because of the very small number of events available, with 
, most of the sources contributing with not more than a single event. 

The analysis of the "clustering" of the observed events can give information about the luminosity of the individual 
sources. After the data of AGASA 0, Q gave hints of possible clustering, this question has received a significant 
. amount of attention 

HIS- 

In this work we want to make a more quantitative and detailed study on how to interpret 

the results on the clustering. 

The nature of the particles in the UHECR remains a central open questions in the field. In this work we will assume 
that most of these particles are protons. This hypothesis is consistent with the existing data (if one takes into account 
the systematic uncertainties on the modeling of hadronic showers) and it is also favored in most theoretical models. 
Here this assumption is also made to allow a detailed quantitative description of particle propagation in intergalactic 
space, which is essential for the problem we are considering. It is possible and reasonably straightforward to generalize 
the discussion to the case of nuclei. 

All cosmic ray sources are unavoidably also sources of neutrinos @, 0] because some fraction of a population 
of relativistic hadrons will necessarily interact with ordinary matter or radiation fields targets inside or near their 
acceleration site creating pions and other weakly decaying particles that can (chain) decay into neutrinos. Viceversa, 
the emission of neutrinos imply the existence of relativistic hadrons, and therefore the acceleration of cosmic rays. The 
catalogues of the CR and neutrino sources are in principle identical. In practice the situation could more complicated, 
because the energy range studied by the CR and neutrino telescopes differ by several orders of magnitude, and the 
relation between the neutrino and cosmic ray emission from a source can vary significantly depending on the structure 
of the source, and it is certainly possible to have bright CR sources that emit a relatively small amount of neutrinos, 
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and viceversa. Nonetheless, it is reasonable to expect that the brighest extragalactic sources of both CR and neutrinos 
will coincide, and it is interesting to discuss in parallel these new astronomies. 

Also in the case of neutrinos one can measure an inclusive flux, that sums the contributions from all sources in the 
universe, while the brightest and most powerful neutrino sources should be identifiable as individual objects. Therefore 
also for neutrinos the "clustering" of the detected events is a useful method of study. An important advantage is that 
for neutrinos linear propagation is certain, but one has to deal with the existence of the foreground of atmospheric 
neutrinos. The theoretical framework needed to discuss the clustering of neutrino events is essentially identical to the 
one needed for protons. A quantitative difference is that for neutrinos extragalactic space is perfectly transparent, and 
the inclusive flux receives most of its contribution from very distant and very faint sources, and the fraction of this 
flux that can be "resolved" in the contribution of identified sources is likely to be small. More that the interpretation 
of the existing data, the goal of this work is the development of some general analysis instruments that can be used 
in the study of future, higher statistics results. 

This work is organized as follows: in the next section we make a preliminary discussion of the so called "Olbers 
Paradox" . The discussion of this celebrated puzzle allows to introduce the key concepts needed to compute the 
signal from the ensemble of all sources in the universe. In section 3 we collect the results on particle propagation in 
intergalactic space that are needed in the following. Section 4 discusses how the flux received from one astrophysical 
source depends on its redshift. Section 5 discusses possible forms of the luminosity function of the proton and neutrino 
sources. Section 6 compute the inclusive particle flux from the combined emission of all sources in the universe. 
Section 7 discusses what part of the inclusive flux can be resolved in the contribution of individual sources. Section 8 
applies these analysis instruments to the interpretation of the recent Auger data. Section 9 discusses extragalactic 
neutrino astronomy. Section 10 gives some conclusions. 



II. THE "KEPLER-OLBERS PARADOX" 

In 1610 Kepler was the first to note that the darkness of the night sky seems to imply that the universe is finite, 
or at least contains a finite number of stars. This point was later rediscussed by several authors, and is commonly 
known in the formulation of Hcnrich Olbers in 1823. The so called "Olbers paradox" is the statement that in an 
infinite, homogeneous and static universe the surface of the celestial sphere should be infinitely bright. A derivation 
of this result is elementary. In euclidean space, the relation between the apparent (I) and absolute (L) luminosity of 
a source placed at a distance R (neglecting photon absorption) is: 
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An homogeneous ensemble of identical sources of density n s produces then an inclusive energy flux per unit of solid 
angle dF/dil that diverges linearly because every spherical layer contained in the interval [R, R + dR] gives an equal 
contribution: 
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Assuming the same homogeneous distribution of identical sources, one can also consider the number of objects in 
the sky that have apparent luminosity in the interval [£,£+ d£]. The distribution dN s /dl can be easily calculated as: 
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a result that is also often stated in its integral form: 
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that has the obvious geometrical meaning of the number of sources contained in a sphere of radius J?^ min = 
y / L/(47r£ m i n ). The infinite brightness of the sky is due to the "exploding" number of very faint sources. In fact 

for l m in - 
diverges as oc 

These elementary results can be also restated in an adimcnsional form, and in terms of particle fluxes instead of 
energy fluxes. We consider identical sources emitting Q particles per unit time, and choose a reference length Rq (here 



the number of observed sources diverges oc £ mi ^ , and correspondingly the integrated energy flux 
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the choice of Rq is completely arbitrary, in the following it will be natural to identify i?o with the Hubble length). 
The particle flux received from a source placed at the distance Rq is: 

$ n = -9— = L ( 5 ) 

Using the notation: R = Rq z and $ = $o % one can rewrite equations (TT|) and ([3]) in the equivalent adimensional 
form: 
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In equations ([6]) and (J7J) we have introduced the two functions F(z) and G(x). The function F(z) = z 2 can be called 
the "dimming function" , and describes how the flux from a source decreases as its distance increases. The function 
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/2 can be called the "flux distribution function" and describes how the number of observable sources 



depends on the size of their flux. 

Summing over all sources, the inclusive flux per unit solid angle (as well as the inclusive energy flux) diverges; 
however it still formally possible to write an expression for the inclusive flux as an integral over the scaled distance z: 



Rq 



(n s Q) 



dz K{z) 



Rq 



(n s Q) 



4tt 'Jo 4tt Jq 



dz 1 — » oo 



(8) 



The function K(z) — 1 can be called the "Kepler-Olbers function" , and gives the contribution of the spherical shell 
between (scaled) distance z and z + dz to the inclusive flux. In general there is an integral relation between K(z) and 
G(x): 



dz K(z) = dx x G(x) 
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(9) 



that is easily understood, noting that both sides of the equation describe the particle flux generated inside a sphere 
around the observer of (scaled) radius z. 

It is obvious that the three functions F(z), G(x) and K(z) are related to each other and, given one of the three, the 
other two can be deduced taking into account the hypothesis of static, euclidean space. For example, if the "dimming 
function" has the general form F(z) the functions G(x) and K(z) can be calculated as: 
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where dV/dz = 4tt Rq z 2 is the volume contained in the spherical shell [z,z + dz}. In fact equations (jTUJ) and (jTT|) 
remain true for a realistic cosmological model, replacing the volume element dV/dz with the correct general relativity 
expression of equation ([55)) . 

The "Kepler-Olbers paradox" is "solved" , and we observe a finite energy flux even if our universe is (to a very 
good approximation) euclidean, (nearly certainly) infinitely large and homogeneous, because the additional condition 
of a static, unchanging universe is false. The calculation of the inclusive flux from the ensemble of all extra-galactic 
sources must obviously take into account the expansion of the universe, and the consequent (redshift) energy loss of 
all particles. In addition the particles can lose energy also for other mechanism that must be taken into account. 
With the inclusion of the cosmological and energy loss effects the functions F(z), G(x) and K(z) take more complex 
forms that are different the simple result: 

F(z) = z-\ G{x ) = ?L11 K(z) = 1 (12) 
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that is valid for a static, euclidean space universe, and in the absence of energy losses. The form of the functions 
F(z), G(x) and K(z) in a realistic cosmological model depend on the particle type, the particle energy and encode 
also the shape of the emission spectrum (and clearly also the parameters that control the expansion of the universe), 
however the solution (|12p will always emerge as the universal asymptotic form for small redshift (z —> 0) and the 
largest fluxes (x — ► oo). The modifications have the consequence that the inclusive particle flux, obtained integrating 
over all possible sources, remains finite. The relation (|9|) remains valid also in the limit of z — > oo (and corresponding 
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The result of the integration in (|13[) is finite and equal to £*■ The physical meaning of can be understood 
qualitatively as the equivalent linear size (in units of the Hubble scale length Rq) of the universe that contributes to 
the flux $, or equivalently as the effective number of Hubble times (to = Ro/c) that contribute to the flux injection. 

III. PARTICLE PROPAGATION 

A. Energy Losses 

In this work we will describe the energy losses of neutrinos and protons propagating in extra-galactic space as 
a continuous deterministic process described by the function dE/dt(E,z) that depends on the particle energy and 
on the epoch z. The assumption of continuous energy loss is essentially exact for neutrinos, that lose energy only 
because of the redshift effect, but is also a good approximation for protons. Integrating the energy loss over time 
it is possible to compute the evolution with time (or redshift) of the energy of a particle constructing the function 
-^evolution (Ei, Zi,Zf) that returns the energy at redshift Zf of a particle with energy Ei at redshift zf. 

Ef = devolution (Ei, Zi, Zf) (14) 

Neutrinos, to a very good approximation, only lose energy for the redshift effect, and their energy evolution is: 

It is particularly important to know the energy at a past epoch z of a particle that is observed "now" with energy E. 
This is a particular case of the general relation (|14p and can be indicated (propagating backward in time) as: 

Eg(E, z) -^evolution 

(E,0,z) (16) 

For neutrinos one obviously has: 

E g (E,z) =E (1 + 2) (17) 

The propagation of protons in intergalactic space has been discussed in several works, after the pioneering work 
of Greisen, Zatsepin and Kuzmin [8j and is now well understood. Together with the redshift losses one has to 
include losses due to interactions with the photons of the cosmic microwave background radiation (CMBR). The two 
important processes are pair production interactions (jr/ — > pe + e~) and inelastic hadronic interactions (with the 
threshold processes p'y — * pir and p'y — ► nir + ). The results of a numerical calculation of E g (E,z) for protons are 
shown in fig. [1] The different curves correspond to the "backward in time" evolution of protons detected "now" with 
E between 10 17 eV and 10 22 eV. Note that the time evolution of the proton energy is completely independent from the 
structure and intensity of the magnetic fields, because these fields can only bend the trajectory of charged particles, 
and the target radiation field is isotropic. Obviously, if the particles do not travel in straight lines, the redshift z has 
no simple relation with the distance of the particle source, but only describes the time of the particle emission. 

Inspection of figure Q] shows that for each observed energy E there is a critical redshift where the energy evolution 
E(z) "explodes", starting to grow faster than an exponential. At this critical redshift the universe, filled now with a 
hotter and denser radiation field, has become opaque to the propagation of the protons. 

The form of the redshift evolution of the proton energy suggests the introduction of the concept of a redshift 
"horizon" . The implicit equation: 



Eg(E, Zh) — £max 
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(with -Emax the maximum energy with which protons are produced in their sources) can be used as a definition of 
the redshift horizon z^E, E max ). An illustration of this definition of the redshift horizon is shown in fig. [51 The 
redshift horizon becomes progressively smaller when the energy E grows. When E becomes a few times 10 19 eV the 
"shrinking" of the horizon with increasing E accelerates sharply. At this "GZK" energy the protons start interacting 
with the photons of the background radiation in the present universe, and not only with the higher energy photons 
of a hotter past. Above this critical "GZK energy" the definition of the horizon with equation (fl8l) strongly depends 
on the "maximum acceleration energy" (that is not a priori a well defined concept) and therefore loses most of its 
usefulness. A better definition of the proton horizon will be given later in equation (|65p . 



B. Magnetic Deviation 

Protons are charged particles, and therefore their trajectories are bent by magnetic fields. Astronomy with protons 
is therefore possible only in a limited region of source distance and particle energy so that the observed deviations 
are acceptably small. The size of this region is determined by the strength and structure the astrophysical magnetic 
fields. 

Neglecting the deviations due to to magnetic fields in the vicinity of the source and in the Milky Way halo, the 
deviation of a charged particle traveling in intcrgalactic space for a distance d can be described (in the limit of small 
deviation) with the expression 

Wl*» - (^f) 2 = (£) (^) 2 = [ZeBf X B (A) (19) 

where rj, — E/(ZeB) is the Larmor radius, B is the typical value of the extra-galactic magnetic field and is 
its coherence length. In (|19[) the magnetic deviation is modeled as the incoherent sum of contributions received in 
"magnetic domains" , where the field has different orientation. In this simple model, the angular deviation accumulated 
by a proton of energy E from a source at distance d scales oc y/d/E, with an absolute value controled by the combination 

B^/Xb~. 

If we require an image with an angular resolution better that 69, for a given energy E we must limit our observations 
to distances smaller than an "imaging radius" : 

E 2 69 2 

Ri^n S (E,S9)= {ZeB)2 ^ (20) 
Alternatively to produce an image with resolution 66 of a source at a distance d one must select a threshold energy 

-^imaging 

£imaging(d, 66) = ^| (Z e B) (21) 

OU 

The combination B^/Xb is therefore the crucial parameter that controls the perspectives of extragalactic cosmic 
ray astronomy. The best (and perhaps today the only) way to estimate B 2 X b and determine where CR astronomy is 
possible is in fact the observation of UHECR and the use of equation (fT5)) . At present however the estimate of B 2 Xb 
remains very poorly known. 

The recent results of the Auger collaboration BUI indicates that the directions of the highest energy events (E > 
5.6 x 10 19 eV) are correlated with the positions of close AGN's, selected with a redhift cut z < 0.018 that corresponds 
to d ~ 75 Mpc, using a correlation cone of opening angle ip ~ 3.1°. The simplest ("nominal") interpretation of these 
results is that the highest energy CR are in fact produced in these AGN (or in sources very strongly correlated in 
position) and arrive with a magnetic deviation of less than 3.1°. This is a remarkable result that, attributing all of the 
observed deviation to the intergalactic propagation (that is neglecting the contributions of the source envelope and 
of the Milky Way) , allows to set an upper limit on the combination B 2 Xb that describes the extragalactic magnetic 
field: 

Z 2 B 2 X B < 0.15 (nGauss) 2 Mpc (22) 
Assuming that the particles are protons, this implies the imaging radius: 

E \ ( $$imagc 



Paging = 210 ( ) ( ) Mpc (23) 
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This result, if confirmed, indicates that CR astronomy should be able to give important results in the very near future. 

The "nominal interpretation" of the Auger results, suffers from several difficulties both internally (see for example 
the criticism in fill ) and comparing with the observations of other detectors (for the Yakutsk data see [l2j for the 
HiRcs data see [la])- The consensus is that the Auger results do indicate the existence of an anisotropy of the arrival 
directions of the highest energy cosmic rays, but the origin of the anisotropy has not been established, the significance 
of the 3.1 degrees of the opening of the correlation cone is therefore ambiguous and the estimates (f22|) and (|23|) only 
tentative. 

A remarkable fact of the Auger results is that two of the highest energy events are in correlation (using the 3.1° 
cone) with the nearby (d ~ 3.5 pc) AGN Cen A. Taking into account the fact that Cen A has been often discussed 
as one of the most likely sources of UHECR this is a striking result, that suggests that the first extragalactic CR 
source has in fact already been imaged. There are also indications that Cen A could be the source of a larger number 
of events with a larger angular spread (see for example the remarks in [II], [l3], UH). For example a 20 degrees cone 
around Cen A contains 9 events. Such a large contribution from Cen A could account for a significant fraction of 
the observed anisotropy effect. These observations could however also a dramatic impact on the estimate of the 
extragalactic magnetic field, because they suggest a larger deviation (59 ~ 10°) for a shorter distance. The estimate 
for the combination B 2 Xb becomes therefore larger by a factor of order 200: 

Z 2 B 2 X B ~ 40 (nGauss) 2 Mpc (24) 

with a correspondingly smaller imaging radius, and significantly restricting the potential of CR astronomy for more 
distant sources. Additional data should allow to clarify this ambiguous situation. The discrimination between the 
two estimates (|22p and (IM1) for the parameter that describes the extragalactic magnetic field is obviously crucial for 
the future of the field. 



IV. FLUX FROM AN INDIVIDUAL SOURCE 



Let us consider a source that emits isotropically particles with the energy spectrum described by the function q(E). 
The emitted particles travel, at least in good approximation, along geodetic trajectories losing energy continuously. 
Propagating backward in time, a particle that is observed "now" with energy E, at the time that correspond to 
redshift z had energy E g (E, z) . 

The flux received from such a source placed at the distance that corresponds to redshift z can be written as: 



Ml 

where: 



Ro = ^- (26) 
is the Hubble length, and the adimensional function F^,(z, E) is given by: 

F (7 F\ — 1 l[E a {E lZ )] dE g (E,z) 

E) r 2 (z)(l + z) tfE) dE— ■ (27) 

In equation ([27)) E g (E, z) is the emission energy (at the source) of the particles observed with energy E, and r(z) is 
the (adimensional) comoving coordinate that corresponds to the redshift z: 

r (*)=l'wi\ (28) 

Jo mz 1 ) 

with TL(z) = H(z)/Hq the (scaled) Hubble parameter at the epoch of redshift z: 



u(z) = y/n m (1 + zf + n A + (1 - n m - ao (1 + zf. (29) 

The expansion of the universe, and the redshift dependence of the Hubble constant are determined by the parameters 
Cl m and 

The "dimming function" F^,(z,E) describes how the flux of a source weakens, and (in general) has its energy 
spectrum distorted when the redshift z increases. The function is the generalization of Gauss law (l/z 2 behaviour) that 
is valid for static, euclidean space in the absence of energy loss. For small z cosmological and energy loss effects become 
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negligible and the function F l j > (z,E) takes the simple asymptotic form F c f ) (z,E) — > z~ 2 . The form of the function 
F^^z, E) is determined by the cosmological parameters (f2 m and f^) an d from the shape of the emission spectrum. In 
principle an unambiguous notation should include these additional dependences: F<f,(z,E) — > F<j,(z, E; {q}, O m , Oa)- 
For simplicity, in this work we will leave the dependence on the cosmological parameters and the shape of the emission 
function implicit. 

A relation analogous to (|25|) can also be written for the integral spectrum: 

$(£ min , z) = Q< f^ F*(z, E min ) (30) 

47T JXq 

where Q(E m i n ) is the integral source emission: 

/>oo 

Q(E min ) = / dE q(E) . (31) 



The function Fg>(z,E) (for the integral flux) can be obtained from F c j ) (z,E) (for the differential flux) with the inte- 
gration: 

F*(z 1 E aia )= 1 / g(£)F (z,£) (32) 

For z — > 0, the asymptotic behaviour of F$(z) takes again the simple form z~ 2 . 

For neutrinos, and more in general in all cases when the only significant source of energy loss is the redshift, the 
energy evolution of a particle has the simple form E g (E, z) = E (1 + z) the general expression for F^^z, E) given in 
equation (|27|) can be simplified, obtaining: 

1 q[E(l + z)\) 

F f {Z ' E) -r^zj q(E) (33) 
As a consistency check one can observe that a source of such "non interacting" particles of absolute luminosity 

L= dE E q{E) (34) 
Jo 

placed at a distance z, has an apparent luminosity: 



/ dE E <j>{E) = [ dEE (— L, q[E(l + z)] 
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This expression has the well known form £ — L/(47tg? 2 ), where c?l is the "luminosity distance": 

d L (z) = Ro r(z) (1 + z) . (36) 

The "dimming functions" for neutrinos take a very simple form if one makes the assumption that the emission 
spectrum is un unbroken power law of slope a. In this case the functions Fj,(z) and F$(z) become independent from 
energy, and equal to each other: 

F 4 ,(z)=F*(z)= {1 + * ) r (37) 
r z (z) 

In case of a power law emission of neutrinos, the spectrum mantains its (featureless) shape even after propagation 
and energy loss. The function F$(z) for neutrinos is shown in figure [3] The function depends on the value of the 
cosmological parameters Q m and Oa, and of the value of the slope a. One can see that F$(z) coincides with the form 
z~ 2 for z < 0.1, and then begins to fall more rapidly suppressing the flux of very distant sources. 

The "dimming function" for proton sources depends very strongly on the proton energy. Some examples of this 
behaviour (for the integral flux) are shown in fig. [4] where the redshift dependence of the function F<f,(z, E m [ n ) is 
shown for several value of the threshold energy E min , assuming again a power law emission spectrum. One can clearly 
see that a proton source becomes rapidly unobservable when it becomes more distant than a rather sharply defined 
"horizon" that shrinks rapidly with energy for E m i n ~ 10 19 to 10 20 eV, and is much smaller that the Hubble length. 
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V. LUMINOSITY FUNCTION 



In this work we will make the assumption that all astrophysical sources have emission of the same shape and differ 
only for their absolute luminosity. The idea behind this assumption is that the emission mechanism is likely to be 
universal, independent from the properties of the individual sources. An example of this situation is the acceleration 
of cosmic rays near the blast waves of SuperNova (SN) explosions, where it is believed that particles are accelerated 
with a power law spectrum oc E~ a with a slope a ~ 2 that is essentially independent from the properties of an 
individual SN (for example the mass and initial velocity of the ejected material), but it is determined only by the 
(universal) compression factor of the interstellar gas after the passage of a strong shock front. 

This assumption of a "universal" spectral shape is of course an hypothesis that could very well be incorrect. This 
possibility is in fact likely to the UHE proton sources, since it is in fact natural to expect that different sources have 
different high energy cutoffs, and that these cutoffs are in the energy range where the spectra are measured. After 
these words of warning, we will anyway make the assumption of a universal shape as it clearly significantly simplifies 
the discussion, and represents a natural first order approximation. 

The luminosity function dn s / dL(L, z) gives the number density (in comoving volume) of sources of luminosity L at 
the epoch of redshift z. The power density C{z) of the sources at redshift z is given integrating over all luminosities: 

C(z) = J dLL^(L,z) (38) 

In this work we will discuss emission spectra that are well represented (at least in the region of interest of the 
measurements) as simple power laws of constant slope. In this situation however, the source bolometric luminosity is 
not unambiguously defined because the integral 

oo 

dE E E~ a 



o 

is always divergent (at the upper limit for a < 2, and at the lower limit for a > 2). In other words, an unbroken power 
law emission is impossible, and to estimate the bolometric luminosity of a source it is necessary to specify the range of 
validity of the power law form, and the behaviour of the spectrum outside this region. This requires the introduction 
of additional parameters that are not easily constrained, because one observes only a limited range of the spectrum. 
What is really needed for our purposes is to determine the absolute normalization of the power law spectrum in the 
energy region of interest. This can be done considering the luminosity in a fixed, relevant energy range. In this work 
we will use the convention to use the luminosity (or luminosity density) integrated over one energy decade starting 
from a minimum energy E m i n that is relevant for the observations 

L = £dccade(£min) - / dE q(E) (39) 



(in the case of a = 2: the luminosity per energy decade is constant and independent from £mm) : To keep the notation 
simple in the following the subscript "decade" will be omitted, and the dependence -Emm (in most cases) will be left 
implicit. The integral emission of particles above the energy i^mm is related to this definition of the luminosity by the 
relation: 

Q(E min ,L) = -^L (40) 

-^min 

where r\ is an adimensional factor that depends on the shape of the spectrum. For a simple power law of slope a one 
has: 

, = ^ = ^(l-10-(-»)) . (41) 

For a — > 2 the limit of equation (|41| is 772 = 1/lnlO ~ 0.434. The differential and integral spectra for an arbitrary 
value of the energy (and assuming the validity of the power law behaviour): 



q(E,L) = 



L rj a (a — 1) 



E. 
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Q(E,L) = 
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(42) 
(43) 
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The simplest hypothesis for the luminosity function is to assume that all sources are identical with luminosity L§: 

^ = S[L - Lo] (44) 

This is nearly certainly incorrect. As an example, the luminosity function of AGN in the hard X-ray band 
([10,20] KeV) has been fitted by Ueda [l6[ with the functional form: 



dn s n ( L c 

~dL [L) -T c It 



r \ 0.86 / \ 2.23 



-1 



(45) 



with L c ~ io 43 ' 94 erg/s. This shape provides a good fit for all luminosities L x £ [10 415 , 10 46 5 ] erg/s, in a broad 
range that extends for 5 orders of magnitude. In the following we will use expression (|45p (shown in figure O as a 
"template" (with two free parameters no and L c , that in principle can also be redshift dependent) to build examples 
of the luminosity functions for the (yet undiscovered) proton and neutrino extragalactic sources. 

Note that the integration of the expression (|45|) for the luminosity function diverges at the lower limit (corresponding 
to a divergent number of faint sources) while the power density remains finite. This fact illustrates the point that 
the concept of "total number of sources" should be treated with caution, because the existence of a large number of 
weak sources is difficult to establish (or refute). In the following we will use the form of the luminosity distribution 
(|45ll assuming also a low luminosity sharp cutoff: 

L > 3 x 10~ 3 L c (46) 

In this work we will use as the "characteristic" value of the source luminosity the median luminosity Lq of the 
sources at redshift z = 0, defined as the solution of the equation 



f Lo dn s , , f°° dn s 

I ~dL = J u it 



(1,0) . (47) 



Half of the power density is due to sources with L < L , and the other half to sources with L > L . For the distribution 
(l45|) with the cutoff (gBJ) one has L = 0.116 L c . 
Without loss of generality, it is possible to write the luminosity function in the form: 



dn S( £(Q) ( L 



^r(L,z) = ^f L (^,z) (48) 



with £(0) the power density of the ensemble of all sources at z = 0, and the adimensional function /l(j/, z) satisfies 
the condition: 

'dyyf L (y,z) = ^ . (49) 

The redshift dependences of the power density of the extragalactic cosmic ray and neutrino sources is of course not 
known, however it is likely that they have behaviours that are similar to each other, and that are also similar to the 
evolution with cosmic time of the star formation rate, or of the AGN activity. It is therefore expected that the power 
density C(z) will grow with increasing z, before falling of at very large redshifts that correspond to early times before 
the formation of the first sources. As an example of possible redshift evolution, in figure [6] we show the evolution of 
the AGN luminosity in the hard X-ray band ([10,20] KeV) as fitted by Ueda [II]. 



VI. INCLUSIVE FLUX 



Summing over all sources, it is possible to define the inclusive particle "emissivity" qi n d(E, z) that gives the number 
of of particles of energy E emitted per unit time and unit of comoving volume at the epoch of redshift z. In the 
simplified model where all sources have emission of the same spectral shape, and the "median source" (of luminosity 
Lo) has emission q(E), the emissivity takes the factorized form: 

qincl (E,z)=q(E) fdL^(L,z)^ = ^lq(E) (50) 
J dL L Q L 
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From the emissity is is possible to obtain the space averaged particle density at the present epoch (n(E)) as: 
(n(E)) = J dE Q J ° dt q hlcl (E , t) S[E - E cvol (E , t, to)] 



q inc i(E,0) f°° 1 qiDC i[E g {E,z),z] dE g (E,z) 



H J H(z){l + z) q in cl(E,0) dE 

The first line in (|5ip simply says that particles observed now with energy E have been produced in the past at time t 
with initial energy E ; the integration over all possible times and energy is limited by a delta function because we are 
assuming that the evolution with time of the particles energy is deterministic, and therefore for each emission time 
t, only particles with a well defined Eo(E,t) can contribute to the density of particles with energy E. The result of 
equation (|5ip can then be obtained transforming the integral over time in one over redshift using the Jacobian factor 

d l = 1 _ 1 1 f52) 

dz H(z) (1 + z) Ho H(z) (1 + 2) v ' 

and performing the integration over the emission energy Eq with use of the delta function. 

If we make the assumption that all sources have identical spectral shape, and therefore that the emissivity is 
factorized with the form (|50[) the we can recast equation (|51[) in the form: 



(„(£)) = g^i(£,0) f 00 dz 1 q[E g (E,z)} C(z) dE g (E 7 z) 



Ho Jo H(z) (1 + z) q{E) £(0) dE 
q ind{ E,0) r dz ^l F ^ E) m (53 ) 



Ho Jo n(z) vy ' ; £(0) 

where we also have used the definition of F${z, E) in equation ([27]) 

It is important to stress the fact that this result describes the space averaged particle density. The particle density 
"here" (meaning near the Earth, or in the disk of the Milky Way) is of course a directly measurable quantity, and does 
not necessarily coincide with the space-averaged one. The "danger" is that intergalactic (or source envelope) magnetic 
fields can be sufficiently strong so that the emitted particles can only diffuse away slowly from their sources. The 
extragalactic cosmic rays could therefore form slowly expanding bubbles around their sources, and be homogeneous 
only after averaging over a distance scale larger than the average source separation. Since the diffusion is a function of 
rigidity this situation could result in a variety of different interesting possibilities, depending on the energy dependence 
of the diffusion coefficients and the position of the observer with respect to the sources. If the particles during their 
lifetime have the time to propagate for a linear distance comparable to the source separation then the extragalactic 
cosmic ray gas becomes homogeneous, these complications can be safely neglected, and we can identify the density 
of the observable extragalactic cosmic rays with the space averaged one given in equation (|53p . This same condition 
on CR propagation also garantees the isotropy of the particles, and the observable inclusive particle flux (f)i nc i(E) (in 
units (cm 2 s sr) _1 ) can be estimated as: 

0inci(£) = ^ (n(E)) (54) 

4 7T 

It is useful to introduce notations for the integral (and the integrand) over the redshift z in equation (|53|) : 

dz St4 F *( E > z ) 7^ = / dz K ^ z > E ) = ^(E) (55) 



H(z) vy ' ; £(0) 
With these definitions the inclusive flux can be written in the form: 

<f>md(E) = ^ q- mc i(E, 0) Q(E) (56) 

The function K ( p(z,E) in equation (|55[) has the meaning of the "redshift response", and gives the contribution of 
the time interval that corresponds redshift in [z, z + dz] to the particle flux at energy E. For small z the function 
K^(z) takes the the asymptotic value K^, — > 1. The quantity (<j>(E) can be understood qualitatively as the number of 
Hubble times that "count" in producing the observed inclusive flux. In the case of a "non-interactive" particle like 
the neutrino C,<t>{E) is clearly of order unity because the age of our universe is in fact of order of one Hubble time, 
In case of strong source evolution, with a much larger activity in the past, the factor C<j> can be as large as 2-3. For 
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ultra-high energy protons the value of C<t>(E) reflects (in units of the Hubble time Hq 1 ) the characteristic time for 
energy loss of the particles and < 1, varying also rapidly with energy. 

The appearance in equation (|53p of the function F^(z, E) (that is dimming function for a source calculated assuming 
linear propagation for the emitted particles) may appear surprising. It can therefore be instructive to give a second 
derivation of the equation (|53|) assuming that the particles propagate along straight lines. In this case, if the particles 
move with a velocity close to the speed of light, and the distance between sources is small with respect to the Hubble 
length c/Hq the condition of homogeneity of the emitted particle density is clearly satisfied. For particles propagating 
at the speed of light along straight lines the redshift z of the emission of a particle is not only in a one-to-one relation 
with cosmic time, but also with the source distance. The inclusive flux can therefore be computed summing the 
contributions of the sources contained in different layers of comoving volume. The inclusive flux is then obtained with 
the double integration: 

<t>inci(E) =-L J^dz^ JdL ^(L,z) tf>{E,L,z) (57) 
where dV/dz is the comoving volume contained in the layer [z, z + dz] 

dn s /dL{L,z) is the number density of sources with luminosity L at the epoch z, and (j)(E,L,z) is the flux from a 
source of luminosity L placed at distance z. For a universal emission shape q(E) (normalized to have the median 
luminosity Lq) one has: 

with F^(z,E) the dimming function defined in (f2"T)l . It is then simple to see that one reobtains expression (|53[) . 
confirming the consistency of the treatment. 

It is also possible to use an analogous notation for the inclusive integral flux introducing the functions K$(z,E) 
and C*(^) : 

3>incl(-Emin) = ~T~ Qincl(^min, 0) C$(-^min) = i. ~ ~ C*(-^min) (60) 

4 7T 4 7T E m \ n 



vM£\ = / </•• lu.f..E, )=/ dz ^| F${z,E min ) (61) 



with 



where the dimming function -F$(z, -E m in) for the integral flux was defined in equation (132)) . 

Some examples of the redshift response (or Kepler-Olbers) function for neutrinos emitted with a power law spectrum 
of slope a {q{E) oc E~ a ) are shown in figure [7] In this particular case (no energy loss except for the redshift effects, 
and a featureless injection spectrum) the redshift response is independent from energy and one has: 

*,(,, = *.( 2 )=<!+ip «i> (62) 

As a consequence also the "shape factors" C<f> = C* are equal to each other and independent from energy. This means 
that an injected power law spectrum of neutrinos is not deformed by propagation effects. As numerical examples, for 
a matter dominated cosmology (f2 m = 1, Q\ — 0) and no source evolution one has: 

2 

C*>(a) = C*(«) = 7T, — ( 63 ) 
I + a 

Softer spectra (larger slope a) have a smaller value of C,^. For the concordance cosmology (Cl m = 0.3, fl\ = 0.7) the 
shape factors are larger by approximately 30%, for example for a = 2 and no evolution = 0.5326. The inclusion 
of source evolution can on the other hand increase the value of the shape factor by a significant factor. For example 
using the same redshift dependence of the AGN X-ray luminosity of shown in figure [6] the shape factor for neutrinos 
increases by a factor ~ 4 (for a — 2 one finds C$ — 2.2). 
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The redshift response for protons has a much more complex and interesting structure. Figure [8] shows the redshift 
response function K^(z, E) for values of the proton energy. One can see that for sufficiently low energy E < 10 17 eV, 
the function takes the same form as for neutrinos, since the role of jry interactions becomes negligible. With increasing 
energy the range of redshift that contribute to the flux becomes progressively smaller. Note also how the energy loss 
effects enhace the role of a narrow range of redshifts where the high energy tail of the emission rapidly lose energy. 

Figure [9] shows the energy dependence of the factor ^{E). For energies much smaller that E ~ 10 18 eV the factor 
C<f>(E) reaches an asymptotic value that is in fact equal to the one for neutrino, while at larger energies the effect of 
C<j>(E) is to distort the power law injection spectrum. The factor (<p(E) plays therefore the role of a "shape factor" 
that marks the "imprints" on energy loss on the observable energy spectrum. Berezinsky and his collaborators [HI 
have discussed how this spectral distortion distortion could be the origin of the "ankle" structure observed in the CR 
flux. 

Figure mil shows the Kepler-Olbers function K$(z, E m i n ) for the integral flux for two values of -Emm (-Emm = 
6 x 10 19 eV and 10 20 eV); the proton injection spectrum is a power law with slope a = 2.0 and 2.4. At these very 
high energies only a small volume of the universe is observable, and the size of this observable volume changes very 
rapidly with the threshold energy. Integrating over z one obtains (for a — 2.0) C* = 0.037 for E min = 6 x 10 19 eV 
and = 0.0107 for E m i n = 10 20 eV. For the softer spectra (a = 2.4) shrinks by approximately 10%. 



The calculation of the inclusive flux as an integral over redshift: 

4>(E) = p- q incl (E) / dzK 4) (z,E) = -^q incl {E)^{E) 
4?r J 4tt 

suggests the most useful definition of an energy dependent redshift horizon using the implicit equation: 



(64) 



r 

Jo 



hor 
'f 



dz K <j> (z,E) = f (tiE) (65) 



For example Zq°q(E) is the redshift that contains the emission of 90% of the protons that are observed with energy 
E. Note that this definition of the horizon depends on the assumed shape of the emission spectrum. 

Figure [15] shows the 90% redshift horizon as a function of the proton energy, assuming injection with a power law of 
slope a. In the same figure we show also the energy dependence of the imaging radius (for a 3° resolution) calculated 
for the two estimates of the parameters B\/\b of equations (|22|) and (|24|) . Proton astronomy is going to be limited 
to a region z < zi wr (E) and z < z ma gnctic(i?) — -Rimaging (E)/R . 



VII. THE INTENSITY AND MULTIPLICITY DISTRIBUTIONS 



If the universe is filled with a homogeneous distribution of sources, an observer will see an ensemble of objects with 
a broad distribution of flux intensity, with the faintest (brighest) objects corresponding to the most distant (closest) 
sources. 

In the following we will consider the the integral flux above a fixed threshold energy E m i n . In the following 
expressions the dependence on this minimium energy will be left implicit. If all sources have the same spectral shape, 
the distribution dN s /d$ that gives the number of sources observable with flux in the interval [$, $ + <i$] can be 
written in the form: 



where the adimensional constant N-r u 



bblc- 



iVHubble = 47ri^ ^p- (67) 

is proportional to the number of sources in one Hubble volume, and the flux (f>o : 

, _ Qo _ Lq 7] 

*° 4nRl 4^i? 2 £ min 108 > 

is the flux from a "typical source" (chosen here as the source with the median luminosity Lq) placed at the Hubble 
distance and calculated neglecting all cosmological and energy loss effects. The general expression for G$(a;) is: 

r°° r 2 (z) 

G»(s)=/ dz-^4- / dyf L (y,z) S[x- y F*{z)] . (69) 
Jo n \ z ) Jo 
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In this expression r(z) and TC(z) depend on the cosmological parameters, F$(z) is the "dimming function" (defined 
in equation (|3"0|) ) and z) is the (scaled) luminosity functions (see equation (05])). The function G$(x) encodes 
the shape of the emission spectrum, the properties of particle propagation and the luminosity function and cosmic 
evolution of the sources. 

As a simple check of this expression one can compute the inclusive particle flux is obtained integrating over all 
fluxes $ the distribution: 

1 f°° dN f°° 

*inci = -r- I d§ $ — i = iVHubbic $o / dx x G 9 (x) (70) 

This expression is identical to the result (|60p because the integral over x x G$(x) is equal to £$. In fact, using the 
definition of G$(x) of equation ([55)) one finds: 

l-oo r 2/ z N / roc 

dx x G$(x) = J dz j—r F&(z) I J dyyf L (y,z) 

™ d -~m F ° {z) m^° (71) 

To be in closer contact with the observations, we will consider an explicit model of detector, described by the 
effective area A(fl) and a finite observation time t and study the number of sources N^ detected with exactly k (an 
integer quantity) events. The multiplicity distribution Nk can be obtained starting from the function dN s /dm, where 
m (a real, continuosly varying quantity) is the expected signal from a source. The quantity dN s /dm gives the number 
of sources with expected signal in the interval [m,m + dm]. The probability to detect k events from a source with 
expected signal m = (k) can be obtained using Poissonian statistics. The expected number of events with multiplicity 
k is therefore: 

f°° ( m m k \ dN s , 
the number of multiplets with multiplicity k > 2 is: 

f°° dN 
JVmuiti = N >2 - / dm [1 - e~ m (1 + m)} —± (73) 
Jo dm 

and the total number of detected sources is: 

f°° dN 
N somccs = N>i= dm [I- e- m ] — ± (74) 
Jo dm 

For the calculation of the function dN s /dm we can observe that the expected signal m from a source of luminosity 
L, redshift z and celestial coordinates Q is: 



m(L, z, Q) = $(L, z) A(n) t = * i = - 4% F*(«) * (75) 

^0 4 7T itg £/ m ; n 



The complication here is that in general the detector exposure depends on the celestial coordinates. A simple modeling 
of the average effective area of the Auger detector and of neutrino telescopes as a function of celestial declination is 
discussed in appendix [AJ We will use the notation: 

A(Q) = (A) a(Q) = — o(fi) (76) 

4 7T 

where (A) is the detector average effective area, Aq is the angle integrated effective area (in units cm 2 sr) and <z(f2) 
is a dimensionless function that satisfies the normalization condition: 

4 7T 

The distribution dN s /dm can then be written as: 

m = aw Gm f -\ (78) 

dm mo V m o / 
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with -/Vfjubbic defined in equation (|67[) and: 



m = 4> (A) t = Lo 2 l (A) t (79) 

4 7T K fc m in 



The "scaling structure" of equations (|66|) and (|T8[) can be easily deduced with simple dimensional analysis observing 
also that the number of sources must be proportional to the quantity C(0)/Lq that is itself proportional to the source 
density, and that the signal from any source is linear in the observation time t, the detector size (proportional to (A)), 
and the source luminosity L. 

The function G m (x) can be calculated as the convolution: 

g -w=/Shp) g *(h|)) < 8 »» 



Another way to express the same result is: 

r 2 {z) 



f dQ f°° 
GM = J - J dz 



H{z) 

The function G m (x) satisfies the "sum rule" : 



dy f L (y,z) S[x-ya(Q) F^(z)] (81) 



f°° ( f d£l \ f f°° \ f°° 

/ dx x G rn (x) — ( / — a(Q) ) I / dx x G®(x) ) — / dx x G$(x) = 
Jo \J 4tt J \J J J 

As a consistency check one can compute the total number of detected events as: 



(82) 



Events = / dm m — — = iV Hu bbio m Q / dx x G m (x) 
am Jq 

RoC^U c \ Ant = $. nclAnt (83) 

4 TT E min J 

The result is simply the inclusive flux (already calculated in equation ([50)1 ) times the total detector exposure An t, as 
it must be. Note that the total number of detected events does not depend on the luminosity function of the sources, 
but only on the (space averaged) emissivity, and also that the result is independent from the geometry or "field of 
view" of the detector. 

The shape of the functions G$(x) or G m (x) in general is complicated, since it encodes the properties of particle 
propagation, the source luminosity function and, in the case of G m (x), also the detector geometry, however the limit 
of large x (that is the largest fluxes or signals) the functions have simple asymptotic behaviours: 

G*(x) — ^ (y 3 / 2 ) (84) 



and 



where the quantity 



depends on the luminosity function, and 



G m (x) — (y 3 / 2 ) (a 3 / 2 ) (85) 



(y 3 / 2 ) = / dyy^ 2 f L (y,0) (86) 



(a3/2) = i o8/9(n) (87) 



depends on the detector location and geometry. These physically intuitive results can be easily obtained computing 
the euclidean/static limits of equations ([69]) and (f8Tj) with three steps: (i) use the asymptotic (small z expression) 
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r 2 (z)/H.(z) — > z 2 ; (ii) approximate the function F$(z) —> z 2 ; and (iii) neglect the redshift dependence of the scaled 
luminosity function Jl{v, z) — > The integral over redshift can then be performed using the delta function, and: 



dz 



z 2 



2 a 



,3/2 

V ^' V2 (88) 



The m -5 / 2 ($ -5 / 2 ) behaviour of the distribution of expected signal (flux) remains valid only for sufficiently large 
signals (fluxes), then the distribution stops growing so quickly for decreasing m, and the corrected distribution is 
integrable down to m — > 0. 

The number of sources with and expected signal larger than the minimum value fh can be calculated integrating 
the distribution dN s /dm: in the interval m > fh. If fh is sufficiently large one can use the asymptotic form (|85|) with 
the result: 

This result has a very transparent meaning if one considers the limit of isotropic acceptance and identical source 
luminosity; then the factors (a 3 / 2 ) and (y 3 ^ 2 ) become unity, the quantity C/Lq is simply the source density, and the 
result takes the form 

N s (> m ) _> n s i£ R 3 m (90) 



3 



where 



m ) \ 4 7r E mln m 

is the radius (in static, euclidean space) of the sphere at wich the signal has value m. Note that the radius R m is 
independent from the Hubble distance Ro, as it must be, since it is calculated neglecting cosmological effects. 

Similarly, the expected number of sources that yield a "multiplet" of multiplicity k > 2 can be calculated using the 
asymptotic form of G m with the result: 

f°° dN _ 

Applets = / dm [1 - e- m (1 + m)] ~ ^ N s (m > 1) (92) 
Jo dm 

where we have used the result: 

dm [1-e— (1 + m)] VL^l = (93) 

The expected number of multiplets can be well approximated as the number of sources with an expected signal larger 
than unity, multiplied by a factor y/ir. The limit of validity of the approximation |92|) will be discussed later. 

Some examples of the function G$ (x) calculated for neutrino sources assuming a power law emission, and all identical 
sources are shown in fig. [TTJ The function G<s>(x) is calculated for two different set of values of the cosmological 
parameters (a matter dominated universe and the concordance cosmology) and two different value of the slope a. 
One can see that for x large the function takes the "universal shape" G$(x) — > x~ 5 ^ 2 /2, while for sufficiently small x 
the curve deviates from the "euclidean form" . The small x (or faint source part) of the function is not universal, but 
depend on the cosmological parameters and the shape of the emission spectrum. 

Figure [T2l shows some examples of the G$(a;) function calculated for protons, integrating the observed flux above 
two different threshold energies E m - m — 6 x 10 19 eV and E mul — 10 20 eV. The calculation is performed for an 
emission spectrum E~ 2 with two assumptions about the sources, either all identical sources, or the broad luminosity 
distribution of equation (|45)l . Again, at large x the function takes its universal form x~ b l 2 /2 for a fixed source 
luminosity or x~ 5 / 2 (y 3 / 2 ) /2 for a luminosity distribution. However for protons this euclidean behaviour stops at a 
larger x, because of the effect of energy loss. In fact for the higher energy threshold, the euclidean behaviour is valid 
in a smaller range of x. 

For a qualitative understanding of the role of the function G m ^(x) one can look at figures [T3l and [T4l that show the 
distribution of expected signal m — $ A t for an isotropic neutrino detector. Neutrinos are emitted by an ensemble of 
identical sources of luminosity Lq emitting with a power law spectrum of slope a — 2; the cosmological parameters 
are VL m = 1 and = 0. The injection power density C is taken to the same value, but different calculations assume 
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different values for the luminosity Lq of the individual sources. Changing Lq leaves the total signal unchanged, but 
generates very different distributions dN /dm. Figure [T4l shows the distribution dN/dm is the form m 2 dN/dm versus 
m with a logarithmic scale for m. In this way the area under the cuve is proportional to inclusive signal and the 
shape of the curve is identical to the shape of the curve x 3 G<s> im (x). It is simple to see that the power density of 
the sources change the absolute normalization of the curve, the luminosity Lq "shifts" the distribution to higher or 
smaller signals, with the overall shape determined by the scaling function G$(x) that is calculable from the particle 
propagation properties, the cosmological parameters, and the shape of the luminosity function. 



VIII. INTERPRETATION OF THE DATA ON UHECR CLUSTERING 



It is possible to use the instruments developed in the previous sections to analyse the results obtained by UHECR 
detectors. As an example we will discuss the results of the Auger detector for the first integrated exposure Aq t = 
9000 (Km 2 yr sr) @, [To| . In this exposure the detector has observed 27 events with energy above E t h — 6 x 10 19 eV. 
From this event rate it is possible to estimate the power density £ p (Eth) of the emission in the energy interval 
[E t h, 10-Eth] making the assumptions that the particles are protons, and using equation ()83|) that can rewritten as: 



= 17.2 



An t 



10 36 erg/s Mpc 3 J \ 9000 Km 2 yr sr 
From the observed event rate one can estimate the power density is: 



C*(-Eth) 
0.037 



(—) 

V 0.434/ 



£p(Eth) 



1.6 x 10 



gth \ / 0-037 
60 EeV J U*(£th) 



0.434 
V 



erg 



s Mpc 3 



(94) 



(95) 



The error on the determination of the power density is dominated by systematic effects associated with uncertainties 
in the determination of the energy scale of the measurement. An error in the energy scale affects obviously the value 
of E t h, m addition, the value of the shape factor C*(£th) is a very steep function of energy. For E ~ E t h this energy 
dependence can be approximated with C&(E) oc _E~ 23 . Therefore one can make estimate the uncertainty on the power 
density as: 



A£ p (£ th ) 

£p{Etb) 



0.20 (stat) 



" (¥) 

V / scale 



(96) 



The shape of the CR energy spectrum is also important. It determines the quantity r\ that connects the integral 
emission and the power density (QmclC-Sth) = £(Eth) ij/E t h), and also the precise value of the shape factor (<s>(E t h)- 
If the spectrum is a power law with a slope that varies in the interval a G [2, 2.6] then 77 (see equation (|4ip ) varies in 
the interval r\ £ [(lnlO) -1 = 0.434, 0.500], with a steeper spectrum requiring less power for the same flux. For the 
same range of slope, the shape factor (<s>(E t h) (for E t h — 60 EeV) varies in the interval € [0.037, 0.032], in this case 
a steeper spectrum has a smaller horizon and therefore requires more power to generate the same flux. The effects of 
changing the slope a go in opposite directions for 77 and and therefore in first approximation cancel each other in 
the estimate of the power density. 

The determination of the shape of the extragalactic flux of CR requires the disentangling of the galactic component. 
A crucial problem is the determination of the transition energy where the two components (galactic and extragalactic) 
are equal. If the transition energy corresponds to the "ankle" at E ~ 10 19 eV, the slope of the spectrum is only poorly 
determined experimentally. Theoretical prejudice has traditionally favored a flat slope a ~ 2. If the extragalactic 
component dominates down to much lower energy as advocated first by Berezinsky and his collaborators [IH , the slope 
a is better determined and is significantly larger than 2. If the evolution of the sources is negligible, one estimates 
a ~ 2.6, if the source evolution has a behaviour similar to the Stellar Formation Rate (or the AGN luminosity 
function) one obtains a ~ 2.4. 

If uncertainties in the spectral shape are not very significant in the estimate of the power density in the energy 
decade above 60 EeV, they have a very dramatic effect if one attempts to estimate a total CR power density integrated 
over all energies. Extrapolating the power low behaviour to the interval {E m i n , E max } one can easily compute: 



£p[E m in, £max] = £p[-Eth, ^ ^th] 



n 



(a -2) 



E„ 



Eth 



Eth 



-Q + 2 



(97) 



In the limit of a = 2 this becomes: 



(98) 
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with a weak, logarithmic dependence on the energy limits. If a is large the low energy extension of the spec- 
trum is problematic because the energy requirement grows rapidy with decreasing -E m ; n , in good approximation 
oc (i?min/£ : tii) _a+2 - As an illustration the extension as an unbroken power law of the extragalactic CR spectrum 
down to E m i n ~ 10 12 eV requires a power approximately 9 (2100) times larger than the estimate of equation (|9"5)) if 
one assumes a slope a = 2 (a = 2.4). 

The clustering of the events detected above the threshold energy = 6 x 10 19 eV can give us information about 
how the power density (|95p is subdivided in individual sources if we assume that the propagation of protons of these 
energies for distances of the order of the corresponding horizon is approximately linear. The general method to 
compute the expected number of clusters with multiplicity Nk has been discussed before, and involves the calculation 
of the distribution of expected signal dN s /dm, and the folding with a Poissonian probability. The results of a 
numerical calculation of iV mu i ti (the number of clusters with multiplity k > 2) and of detected sources as a function 
of the median luminosity Lq is shown in figure 1161 The calculation is performed for the value of the source power 
density that reproduces the inclusive rate, and using two models for the luminosity function: a fixed luminosity, and 
the broad luminosity distribution given in equation (|45[) and shown in figure [5l 

The behaviour of the curves N souvccs (Lq) is simple to understand qualitatively. For small Lq the cosmic ray 
signal is given by the contribution of many faint sources, each contributing a small fraction of an event and one has 
-^sources — -^events, with a negligible probability of observing multiplets. Increasing Lq the number of detected sources 
(for a fixed power density) decreases monotonically as the number density of sources decreases, while the volume 
contained inside the horizon remains fixed. 

The behaviour of N mu iti(Lo) as a function of Lq shows a maximum for Lq = Lq ~ 10 erg/s. For a median 

1/2 

luminosity Lq <C Lq the number of multiplets grows oc L according to equation: 




°.T4( , ),u^r . ) (»" s > o») 



3/2 



10 36 erg/s Mpc 3 / V 1040 er g/ s \9000 Km 2 yr sr 



The behaviour AT mu iti oc \JL~q has a very simple explanation. The radius of the sphere within which a source of 
luminosity Lq is sufficiently bright to generate a multiplet grows approximately tx y/Lo, accordingly the volume that 

3/2 

contains multiplets grows oc L Q ' ; on the other hand, for a fixed luminosity density the number density of the sources 

-i t # 3/2 1 1/2 

is n s oc L , with the combined result iV mu i t i oc L L oc L . This behaviour stops when the radius of detection 
of a source becomes comparable with the "horizon radius" of order (i?oC*)- Additional increase in Lq decreases 
the number of the multiplets, since the detection volume remains approximately constant while the source number 
density decreases. The coefficient of proportionality between -/V mu i t i and \[Lq depends linearly on the factors (a 3 / 2 ) 
and (y 3 / 2 ). The first, defined in equation (|76|) , depends on the detector geometry and event selection criteria. For 
the Auger detector location and the cut on zenith angle 9 < 60° one has (a 3 / 2 ) ~ 1.26. The factor (y 3//2 ) = 1 is unity 
if all sources have the same luminosity; for the broad luminosity distribution PS)) (y 3 ^ 2 ) — 1.64. It is intuitive that 
a broad luminosity distribution (with large (y 3 ^ 2 )) yields a larger number of multiplets than a narrow one thanks to 
the contribution of high luminosity sources in the tail of the distribution. 

In the regime of validity of equation (|99p It is possible to obtain Lq with a simple, closed form expression as a 
function of Advents and iV mu iti solving the system composed of equations and : 

T 1 A A ( NmultiA E th [Rq ($(E th )} 2 1 1 

Lo * 1Un I aw: J v (An t) ww 2 (100) 

Substituting the number of events observed in Auger one finds: 

, 39 (AUiiti) 2 ( E th \ ( C$(£th)V erg 



W.5X10- (^jjfj -f (101) 

Note that in the regime of validity of this equation the number of multiplets grows with time as: Af mu iti oc i 3 / 2 , while 
the number of events is A^ even t s oc t, therefore the combination [(A r mu iti/A'events) 2 is a constant. The behaviour 
Lq oc (N mu \ti / N cvcnts ) 2 loses is validity for Lq > 10 41 erg/s, when sources not too distant from the energy loss horizon 
begin to be visible. 

The question of how the 27 highest energy events detected by Auger can be subdivided into the contribution of 
different sources has not a simple answer, and is intimately related to the problem of the the deviations of the particles 
during their propagation and the size of the extragalactic magnetic fields. 
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The "nominal interpretation" of the Auger events, suggests that particles with E > 60 EeV are deviated by 
approximately 3 degrees, and therefore that the clustering of events should be searched for using a correlation cone 
of approximately this size. If one tentatively accepts this conclusion, defining a cluster as a set of events with angular 
separation less than a 3°, one finds 2 doublets and correspondingly 23 singlets. Increasing the size of the correlation 
cone to 5° one finds one triplet and three doublets, while the 18 remaining events are singlets. These numbers should 
be corrected for the effect of random coincidences that in turn depends on the assumptions about the larger scale 
angular distribution of the flux. Taken at face value these results, together with the the calculations shown in figure [TBI 
suggest an allowed range for the luminosity Lq: 

10 39 ^ < L < 10 41 - 5 HI (102) 

s s 

The lower limit follows from the fact that (at a modest confidence level) some clustering is observed. The upper 
limit follows from the fact that most of the events are actually in singlets, and that the total number of objects that 
contribute to the signal is larger than ~ 20. 

In fact the Auger data suggest that at least one source of UHECR has been established with a very high degree 
of confidence. The position of the very close Active Galactic Nucleus Cen A coincides with the most prominent of 
all tentative clustering structures (the triplet found for a 5° aperture cone). Considering that Cen A, on the basis of 
its properties and closeness to the Earth, has been proposed as one of the most likely sources of UHECR this is a 
striking result. The detection of Cen A is an remarkable result of historic importance. It is the first imaging of an 
astrophysical object with charged particles and, 4 centuries after the invention of the optical telescope and 20 years 
after the detection of SN1987A with neutrinos, marks the birth of a new astronomy. 

The clarification of the size and structure of the excess of events from the Cen A direction is very important to 
evaluate the perspectives of the field of CR astronomy. In fact it is possible that Cen A accounts for a number of 
events larger than 2 or 3. For example a 10° (20°) cone around the direction of Cen A contains 4 (9) events, suggesting 
that the excess from this direction is larger and more extended. It should be stressed that Cen A is remarkably close 
to us with a distance of only 3.5 pc. Since the deviations of a particle during propagation in intergalactic space scales 
as Vd (with d the sorce distance) a signal with an angular size of several degrees from Cen A suggest that astronomy 
from sources at larger distances will be very difficult and require higher energies. The result is also in possible conflict 
with the correlation cone of 3.2 degrees found between the direction of the highest energy cosmic ray and the position 
of AGN with z < 0.017 found in Q since the distance of these AGN is typically 20 times larger, and the expected 
deviation of particles produced in sources is expected to be 4 times larger. 

There are several solution to this problem: the apparent large angular size of the event excess in the direction of 
Cen A could be due to the contributions of other sources in the same portion of the sky; the large angular deviations 
could be the effect of strong magnetic fields in the vicinity of the source and not in intergalactic space (but of course 
this should then be expected for a significant fraction of all the other sources) ; the angular spread of the signal could 
be associated with the size of the emission region from Cen A that has jets extending several degrees away from its 
nucleus. It could however also mean that the intergalactic magnetic fields are stronger than expected, and that cosmic 
ray extragalactic astronomy is only possible only for very close sources or very high energy. 

In any case, the power emitted in cosmic ray from CenA is straightforward to compute. In fact the source is so 
close to us that the cosmological and energy loss corrections are small for E ~ 6 x 10 19 eV. The luminosity in the 
energy decade above E m i n can be estimated as: 

£ — (dy 2 ? <™> 

(where -/VconA is the number of events arriving from Cen A). This result is in reasonable agreement (at least in order 
of magnitude) with the estimates (jlOip and (|102[) . especially allowing for a broad luminosity distribution. 

It is interesting to make a prediction on how the number of detected clusters will scale with additional exposure. 
An illustration of this problem is in fig. [T7] that shows how the number of clusters with multiplicity n > 2 grows 
with exposure for the Auger detector. If the median luminosity of the sources is small the number of clusters grows 
cx (At) 3 / 2 , faster than the linear growth of the inclusive rate. Eventually this growth must slow down and finally 
stop when the exposure is sufficient to observe as clusters of events sources close to the energy loss horizon. The 
asymptotic number of detectable sources is clearly sensitive to the form of the luminosity function, and to the density 
of weak sources with luminosity much smaller than the median luminosity Lq. 
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IX. EXTRAGALACTIC NEUTRINOS 

The existence of the extragalactic cosmic rays is the best reason to predict also a measurable flux of extragalactic 
neutrinos. In fact, the CR sources must unavoidably also generate neutrinos, because a fraction of the relativistic 
hadrons they produce will interact with matter or radiation fields inside or near the source, producing weakly decaying 
particles (pions and kaons) that generate neutrinos. The power density and luminosity function of the neutrino 
extragalactic sources are therefore related to the analogous quantities for cosmic ray emission. The relation between 
the CR and the v emissions is very uncertain and depends on the structure of the sources. The motivation for 
the simultaneous measurement of both cosmic rays and neutrinos is precisely the understanding of the nature and 
structure of their common sources. 

In the following wc want to discuss the observations of neutrinos using the detection of ^-induced muons in "Km 3 " 
type neutrino telescopes like IceCube at the South Pole and Km3Net in the Mediterranean sea. We will use a simplified 
model for the effective area of a neutrino detector that is discussed in appendix [B] It can be useful to think about the 
signal of neutrino induced muons as a measurement of the integral neutrino flux above the threshold £™ ln ~ 1 TeV 
using an effective area that is the geometric area of the detector times an energy averaged neutrino to muon conversion 
efficiency of order (e^-^) ~3x 10~ 6 . The exact value of (e^-^) depends on the shape of the neutrino spectrum, 
the energy threshold for muon detection and the average zenith angle of the source considered. Taking into account 
the small conversion efficiency the effective area of a Km 3 neutrino detector is of order of only few square meters, 
many orders of magnitude smaller that the effective area of a telescope like Auger (with a surface of approximately 
3000 Km 2 ). It is however interesting to note that, in some sense, the sensitivities of detectors like Auger and IceCube 
are comparable. This can be illustrated with the following example, if we have one source that emits protons and 
neutrinos with spectra that are approximately equal in both shape and absolute normalization, and the spectra can 
be reasonably well approximated as power law of slope a ~ 2, the the source (if in the field of view of the detectors) 
will leave event rates in ultra high energy showers, and in neutrino induced muons are of the same order of magnitude. 
This can be easily understood qualitatively observing that energy threshold shower detection is of order ~ 10 20 eV, 
while the neutrino telescopes is sensitive to neutrinos with E v > 10 12 eV. The integral fluxes above these energy 
thresholds are (for a spectrum of slope a ~ 2) in a ratio 10 -8 , but also the effective areas differ by approximately the 
same factor, resulting in approximately the same event rates. 

This discussion is made more quantitative in figure[T8]that show the curves in the plane {z, L} (distance, luminosity) 
that corresponds to one event /year in Auger (for a source places at S ~ —30°) and in a Km 3 telescope of the IceCube 
type. The luminosity must be understood as the luminosity per energy decade, and the source spectrum is a power 
law with slope a — 2. For very close sources, the lines that correspond to a constant event rate have the simple form 
L(z)z 2 ~ const. For these close sources, the neutrino event rate is a factor approximately 5 (3) times smaller than the 
shower rate with a threshold £ = 6x 10 19 eV (E = 10 20 eV). On the other hand the proton source become effectively 
invisible when they are behind an energy loss horizon that is a rapidly varying function of energy, while the entire 
universe is visible with neutrinos. 

It is interesting to perform also for extragalactic neutrinos the same study done for protons, and discuss what 
fraction of the inclusive signal can be resolved as the contribution of individual sources. 

With respect to the study of UHECR there are a several obvious differences. The first one is that neutrinos travel in 
straight lines and therefore the problem of the impact of extragalactic magnetic field is entirely absent. The angular 
size of a cluster of j/-induced muons is controled by the muon-neutrino angle, that is calculable (the angular dimension 
of the signal has a weak but non entirely negligible dependence on the shaope of the neutrino spectrum) . 

A second difference is the presence of the atmospheric neutrino background that constitute an important complica- 
tion. The total number of neutrino events from astrophysical sources can only be obtained with a non trivial analysis 
that requires the estimate and subtraction of the background. 

An additional problem is that much (but not all) of the information about the shape of the neutrino energy spectrum 
is lost. The detected muon energy can only be measured with poor resolution, and this energy is only poorly correlated 
with the parent neutrino energy. 

Finally, neutrinos can travel without interacting across a Hubble distance, and it is therefore expected that the 
fraction of the inclusive signal that can be resolved and associated to individual sources, is going to be much smaller 
than for UHECR that can reach us only from a much smaller volume determined by the energy loss length. 

At present there is only an upper limits on the possible flux of extragalactic neutrinos. The most stringent one has 
been obtained by the Amanda detector at the South Pole [l7| , and is expressed in reference to an isotropic neutrino 
spectrum of form </>„ = K v E~ 2 . The 90% C.L. upper limit on the constant K v (relevant for the combination -\-V^) 
is 



K v „ = E 2 <f>^ < 7.4 x 10~ 8 GeV/(cm 2 s sr) 



(104) 
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The quantity K v can be expressed in terms of the power density of the neutrino sources at the present epoch using: 



3^ 



4tt In 10 



where the factor of 3 takes into account the expected contribution of all other neutrino flavors (see [lj 
discussion). The upper limit on the power density of the extragalactic neutrino sources is: 



(105) 



for a critical 



£?<1.04xl0 37 4^ 



2.2 

c 



erg 
Mpc 3 s 



(106) 



(The estimate ~ 2.2 includes an important redshift dependence of the source power density). It can be useful 
to translate this upper limit on the extragalactic neutrino flux as a maximum number of events generated from 
extragalactic neutrinos. The total number of extragalactic neutrino events can be calculated using equations (1831) . 
and depends on the power density of the neutrino events: 
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(107) 



A calculation of the number of sources and of the number of multiplets as a function of the median luminosity Lq 
is shown in fig. [T9l The calculation is performed for a power density of the neutrino sources £° = 10 37 erg/(s Mpc 3 ) 
just below the existing constraints. For Lq not too large the expected number of multiplets is proportional to \/Lq 
according to the approximation© of equation (|99p that can be recast in numerical form as: 
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(108) 



In equation (|108p the quantity (a 3 / 2 ) has been estimated as y/2. This value is valid for IceCube that observes only 
half of the sky. For a neutrino telescope placed in the mediterranean sea at a latitude of approximately 36° one can 



estimate (a 



3/2\ 



1.11, that is 20% lower. For a neutrino detector placed at the equator one should have (a 



3/2\ 



1. 



The dependence of the expected number of resolve sources on the detector latitude is an elementary geometric effect 
associated with the reduced field of view of a detector close to the poles. A detector that observes only a fraction 
1/n of the celestial sphere (with constant effective area in this portion of the sky) has (a 3 / 2 ) = y/n. The "horizon" 
of such a detector is yfn time deeper that the horizon of a detector that observes the entire sky with the same angle 
integrated exposure, and the number of observed sources scales as the volume contained inside the horizon, that is 
■nT 1 x rt 3 / 2 = yfn larger that for an isotropic detector. Deeper observations of a smaller fraction of the sky yield more 
resolved sources. 

The number of multiplets is also proportional to the total number of events. This relation can be written as: 
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This can also be rewritten as a function of the quantity K v = E 2 § v (E) 
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For completeness, figure [20] shows a prediction of the evolution of the number of clusters of neutrino induced muons 
produced by extragalactic neutrino sources with increasing exposure. The number of multiplets is expected to grow 
approximately oc (At) 3 / 2 even for exposures of hundreds of square kilometer year. This is a consequence of the fact 
that this functional form continues until the detector sensitivity is sufficient to see as muon multiplets sources located 
at approximately a Hubble distance, and this does require very large exposures, unless the luminosity of the source is 
extraordinarily large. 
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X. OUTLOOK 



The key problems about the nature and properties of the extragalactic cosmic rays sources remain open, even if 
some of the answers seems to be tantalizingly close. One of the most urgent problems is the determination of the 
region of source distance and particle energy where astronomy with charged cosmic rays is possible. The results of 
Auger on the correlation of the arrival direction of the highest energy CR with the positions of close AGN suggest 
that astronomy with CR is possible with a threshold of ~ 6 x 10 19 eV, and a volume of radius R ~ 75 Mpc. These are 
very encouraging results, but they need to be confirmed with more statistics. If a significant part of the anisotropy 
effect observed by Auger is due to the contribution of a single source at only 3.5 pc the observation of more distant 
sources could be significanly more problematic. 

Assuming that the linear propagation of the highest energy cosmic rays is a reasonable approximation and that 
most of the particles are protons, the clustering (and lack of clustering) of the highest energy events determinee 
to within a factor of 300 or so the typical luminosity of a extragalactic CR sources (in the energy decade above 
E ~ 6 x 10 19 eV) as Lq £ [10 39 , 10 41 5 ] erg/s. This is not inconsistent with the estimated luminosity of Cen A that 
is of order few times 10 39 erg/s. In the same energy region the power density of the ensemble of the CR sources if of 
order of 1 2 x 10 36 erg/(s Mpc 3 ). 

The neutrino observations, that at present limit the power density of the neutrino sources to a maximum value of 
order £„ < 10 37 erg/(s Mpc 3 ) offer already a quite significant constraint, since the neutrino limit refer to the power 
density per energy decade in at much lower energy, and therefore suggest that the spectrum of the CR emission 
cannot extend to low energy with a too steep spectrum. The inclusive extragalactic neutrino flux is mostly generated 
by distant, very faint sources, and therefore only a small fraction (probably less than 10 -3 ) can be resolved in the 
contribution of individual sources by detectors of the IceCube type. If an extragalactic neutrino component is not too 
far from the existing upper limits, the sensitivity of the Km 3 neutrino telescopes could be insufficient to detect more 
than a handful of extragalactic sources. 

Acknowledgments I would like to acknowledge discussions with Steve Barwick, Maurizio Lusignoli, Andrea Silvestri 
and Todor Stanev. 



APPENDIX A: DETECTOR ACCEPTANCE FOR EAS ARRAYS 

The geometry of an Extensive Air Shower Array (EAS) is sufficiently simple that it is possible to have a reasonably 
accurate description of its acceptance with simple analytic expressions. In this type of detection method, cosmic ray 
showers are analysed if their zenith angle is smaller than a maximum value # ma x, and if the shower axis falls inside a 
predetermined region of area Aq . The effective area for the detection of showers coming from a direction with zenith 
angle 9 is therefore A(8) ~ A n cos 9. The angle integrated area is: 

An = A a 2ir f dcos9 cos6> = A tt [l - (cos6» max ) 2 ] (Al) 

J COS ma .x 

The Auger detector has used the cut # ma x = 60°, and therefore its integrated exposure is approximately: 

An — A (^fj (A2) 

The interesting problem is to compute the detector acceptance as a function of the celestial direction. This can be 
done taking into account the trajectories in the detector system of points with different coordinates on the celestial 
sphere. All points with the same celestial declination 6 make the same trajectory in the sky with zenith angle: 

cos 9(h, 6, (pi a t) = cos 6 cos</?iat cos h + sin 6 suiip\ at (A3) 

(h is the hour angle). Averaging over one sidereal day and taking into account only of the time when a celestial direction 
is sufficiently high above the horizon, one can compute the "average" effective area for points with declination 5 as: 

/>27T TT 

(A(S)) = A (cos 9(6, ^iat)) =A — cos 9(h, 6, cp lat ) 9 [cos 9 - cos max ] (A4) 

Jo 27r 

The averaging over a sidereal day is a good approximation for long periods of data taking. The function (cos 9(6, </?i a t) 
for the latitude of the Auger detector ((p\ a t — 24.8°) is shown in fig. [2TJ Integrating the acceptance over all the 



22 



celestial sphere one must obtain the same result as integrating over the detector angular coordinates, and therefore 
one obtains (for any detector latitude) the relation: 

tt/2 x 

d6 cos 6 (cos 0(5, <pi a t)) = — [1 — (cos6> max ) 2 ] (A5) 

APPENDIX B: NEUTRINO DETECTION 

The most promising technique for the detection of point-like astrophysical neutrinos sources is the observations of 
i/-induced muons. These muons are produced in the charged current interactions of Vy, and V ^ in the matter below 
the detector, and are approximately collinear with the parent neutrino direction, and therefore excellent tools for the 
imaging of neutrino sources. 

The z^-induced muon signal is detectable only eliminating the large background of atmospheric muons due to 
the decay of charged pions and kaons produced in the hadronic showers of cosmic rays interacting in the Earth's 
atmosphere. This can be achieved limiting the neutrino observations to up-going trajectories. With this limitation, 
a neutrino telescope placed at the South Pole can only observe the northern celestial hemisphere. 

The signal of astrophysical neutrinos remains partially hidden by the background of atmospheric neutrinos that 
has the same origin of the muon background (the weak decays of hadronic produced by cosmic ray showers). The 
astrophysical neutrino signal can however be disentangled from this background because of its different (harder) 
energy spectrum. Point sources can be identified as clusters of events with approximately the same direction. 

In this work, following [2C| , we will consider a simplified description of neutrino detection that is sufficient for our 
discussion. We will assume that the neutrino energy can be reasonably well described in the relevant energy range 
(I0 12 < E v < 10 16 eV) as a simple power law of slope a. The flux of muons generated by neutrinos with spectrum 
(j) v (E u ) can be calculated as: 

*m(£;T) = / dE » ^ [Psu%(E„e) Y v ^-{E v ,Ef*) 

+ P^(E V , 6) Y v ^ l+ {E v ,Ef»)} (BI) 

In this expression <j) v {E u ) is the neutrino flux summed over all v types, and the factor 1/6 follows from the assumption 
that flavor oscillations mix the neutrinos, so that one has approximately equal amounts for each type (v e , v e v^, Vy 

v r and v T ). For a critical discussion of this result see for example [l9| . The function Psuiv{E v ) (Psurv(E u )) is the 
probability that a v^, (j7 M ) that arrives at the detector with zenith angle 9, survives the crossing the Earth without 
interacting, and Y v ^ pu - {E v , E^ )m i n ) (Yp (E„, ii^min)) is the muon yield, that is the number of muons with 
E > -E^min generated by a (p^) of energy E v . 

The integrand of equation (|BI[) is shown in fig. 1221 for two examples of the slope a and on the two extreme 
cases (minimum and maximum) of the effects of absorption in the Earth. Inspection of fig. [H] shows that for a 
neutrino spectrum of slope not much larger than a ~ 2 the muon signal is generated by neutrinos in the energy 
range [10 12 , 10 16 ] eV. It is therefore useful, at least for a qualitative understanding to consider the flux the neutrino 
induced muons as proportional to the integral flux of neutrino above a (conventional but roughly appropriate) energy 
threshold of E™ m = 10 12 eV, introducing an energy averaged neutrino to muon conversion efficiency {e v ^y) that can 
be calculated performing the ratio between the muon and the neutrino fluxes: 

" " (B2) 



x & V (E U > 1 TeV) 

The conversion efficiency (e^—,^) depends on (i) the muon energy threshold, (ii) the slope of the neutrino spectrum 
(being higher for flatter spectra), and because of the effects of v absorption in the Earth, (iii) on the neutrino zenith 
angle. This is illustrated in fig. [23]) that shows some examples of (e„^ M ). In the following, we will neglect the zenith 
angle dependence of (e^^), considering that this is a reasonable approximation given the size of the astrophysical 
uncertainties. 

For a detector located at the South Pole, the northern celestial hemisphere remains constantly below the horizon 
and can observed continuosly, while the southern hemisphere in unobservable. The effective area will be simply 
modeled as: 

AM Jjjjg (B3) 
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with Aq ~ 1 Km 2 . The angle integrated effective area is Aq — 2ttAq (e„^^). For first order numerical estimates in 
the following we will use = 3 x 10~ 6 , that is the correct value for a low muon energy threshold and a slope 

a ~ 2.1. For a flatter spectra with a ~ 2.5, the conversion is smaller: (e v ^,^) ~2x 10~ 6 . 
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FIG. 1: Redshift evolution of the proton energy. The different energies correspond to protons observed at z = with energy E n 
10 22-o.i,fc eV (with k = 0,1,. . .,25). 
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FIG. 2: Redshift horizon as a function of the observed proton energy E. The different curves are calculated using equation (If 81 for 
different values of the maximum acceleration energy -B max - 




FIG. 3: Plot of the function F^(z) = F$(z) that gives the redshift dependence of the flux from a v source. The source spectrum is a 
power law of slope a. The different curves are calculated for two different choices of the cosmological parameters (O m and Oa)i an d two 
different slopes (a = 2 and a = 2.5). 
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FIG. 4: Plot of the function F^(2,E m ; n ) (in the form z 2 F$(z) versus z) for protons. The function F$(z) gives the "dimming" of the 
integral proton flux above the minimum energy -E m i n with increasing redshift z. The solid (dashed) curves are calculated for a proton 
spectrum slope with slope a = 2.2 (a = 2.6), and different values of the threshold energy E m i n . 
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FIG. 6: Redshift evolution of the X— ray luminosity density of AGN in the [10,20] KeV band according to the fit of [l6f| ■ 
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FIG. 7: Plot of the redshift response ("Kepler-Olbers") function K ( j > (z) = K$(z) for neutrinos assuming a power law emission of slope 
a. In three of the calculations the neutrino emission is considered constant in time. In a fourth calculation the neutrino emission changes 
in time with the redshift dependence shown in fig. [6] 
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FIG. 8: Plot of the redshift response ("Kepler— Olbers") function K^(z,E) for protons assuming a power law emission of slope a = 2.3. 
The different curves are calculated (in the absence of source evolution) for different values of the observed proton energy. The thick (red) 
curve is the curve calculated for neutrinos (considering only redshift energy losses). 
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FIG. 9: Shape factor Ctb{E) for protons, calculated for a power injection with slope a = 2.0, 2.2, 2.4 and 2.6. The solid lines are calculated 
for a constant value of the source power density, while the dashed lines assume for C p (z) the redshift dependence shown in fig. [6] 
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FIG. 10: Redshift dependence of the Kepler-Olbers functions K$(z, i? m i n ) for 

^min — 6 x 10"^ eV (thick curves) and -E m in — 10 20 eV 
(thin curves). The solid (dashed) curves assume a power law injection spectrum with a slope a = 2.0 (a = 2.4). Integrating over z one 
obtain the values C$ = 0.037 (0.034) for E min = 6 X 10 19 eV and £« = 0.0107 (0.0093) for B min = 10 20 eV. 
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FIG. 11: Plot of the function G r f 1 (x) = G$(x) for neutrinos observations. The neutrino emission has a power law form with slope a, 
and all neutrino sources have identical luminosity. The different curves are calculated for different cosmological parameters, and for two 
different values of the slope a. The dotted curve show the function x _5 / 2 /2. 
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FIG. 12: Plot of the function G$(x) (in the form x 2 G$(x) versus x) for the proton integral flux. The calculation is performed for a power 
law energy spectrum with a = 2.4 and two energy threshold E min = 6 X 10 19 eV and E min = 10 20 eV. The dashed lines are calculated 
assuming that all sources have the same luminosity, while the solid line assume the luminosity function of equation 1145 \ shown also in 
fig. [5] The dotted curves show the functions x _5 / 2 /2 and a; -5 / 2 (y 3 / 2 }/2 
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FIG. 13: Number of sources with an expected signal m (plotted in the form m dN/dm versus m). The distributions are calculated for 
the sam, power density £„, assuming that all sources have the same luminosity L and emit with a a power law spectrum of slope a = 2, 
The different curves differ only because of the value of L. 
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FIG. 14: As in figure H"3l but plotted in the form m 2 dN/dm versus m. The area under each curve is proportional to the total number of 
detected neutrino events. 
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FIG. 15: The thick curves show the redshift horizon as a function of the proton energy. The horizon is defined as the value of z within 
which 90% of the observed flux is emitted. The solid (dashed) line is calculated for a power law emission with slope a = 2.7 (2.0). The 
dotted curves indicate the selection criteria used by the Auger collaboration in their correlation analysis . The redshift that corresponds 
to the Cen A position is also indicated. The curves labeld B\ and B2 correspond to an expected deviation of 3 degrees with the estimates 
of the extragalactic magnetic field of equation 11221 and (1241 . 
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FIG. 16: Expected total number of detected sources (solid lines) and number of sources with multiplicity k > 2 (dashed lines) plotted as 



a function of the median luminosity of the sources. The calculation assumes a threshold energy E n 



: 6 X 10 eV and a total exposure 



of 9000 (Km yr sr) for the Auger detector. The source power density is set to C p ~ 1.6 X 10 erg/(s Mpc^) to reproduce the observed 
event rate. The thin (thick) curves assumed a fixed luminosity (the luminosity function of equation (145 D - The dotted line is <x \/ Lq. 
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FIG. 17: Expected number of proton sources detectable in Auger with multiplicity n > 2 (for E min = 6 X 10 19 eV) plotted as a function 
of detector live time. The power density of the sources is C p = 1.6 X 10 36 erg/(s Mpc 3 ). to reproduce the observed inclusive event rate 
(the dotted line gives the total number of events as a function of time). The different curves arc calculated assuming that all sources have 
the same spectral shape (a power law of slope a = 2.0) and 3 possible values of the median luminosity: (Lo = 10 38 , 10 40 and 10 42 erg/s). 
The dashed line are calculated with a fixed luminosity, while the solid lines are the result using the luminosity function of equation l|45j l 
and shown in fig. \5\ 



£:£i 

u 

CD 



O 

■H 



O/i 

O 



46 



44 



42 



40 



38 



36 
10 



_ 


III 


| 1 1 I Mill 

10 20 


II II 

p 


II 1 1 1 1 1 1 1 1 1 
6#10 19 ' / ^ 




•< 




















w 










u 














1 event/year 








Auger South 







1 1 1 III 


Km v Telescope 
mi 1 



-4 



10 



-3 



10 



-2 



10 



-1 



io 1 - 



Redshift z 



FIG. 18: The lines describe the curves in the {z, L} (luminosity, distance) plane that correspond to a rate of one event/year from a 
source. The luminosity must be understood as luminosity per energy decade. The different curves are for the detection of proton showers 
with threshold energy -E m i n = 60 and 100 EeV in Auger (assuming a source declination <5 ~ — 30o), and for the number of neutrino induced 
muons in IceCube (for a source with <5 > 0). The dotted lines are oc z' 2 . 
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FIG. 19: Expected number of neutrino event clusters with k > 2 plotted as a function of the median luminosity Lq of the neutrino 
sources. The calculation is performed with a neutrino power density £^ = 10 37 erg/(s Mpc 3 ), and assuming no cosmological evolution of 
the sources. The solid lines assume the luminosity function of equation ]45t . the dashed lines a fixed value L = Lq. 
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FIG. 20: Expected number of cxtragalactic neutrino sources detectable in IceCube with a multiplicity of up-going muons n > 2 plotted as 
a function of detector live time. The neutrino sources are all identical and emit with power law spectra with slope a = 2 and a luminosity 
per energy decade Lq. The power density of the ensemble of all neutrino sources is C = 10 37 erg/(s Mpc 3 ). The different curves are 
calculated for different values Lq = 10 38 , 10 40 and 10 42 erg/s. The dashed line give the inclusive number of detected up-going muons 
that are produced by extragalactic neutrinos. For a detector places in the Mediterranean sea at a latitude of 36° the number of sources 
should be reduced by a factor 0.787. 
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FIG. 21: Plot of the detector effective area A(Q) as a function of the celestial declination 5. In the case of the Auger— South detector, the 
effective area is divide by the surface area of the array, the detector latitude is <^i at = —24.8° and the projected area is averaged for the 
time when a source has zenith angle 8 < 60°. For the neutrino telescopes, the detector projected area in the direction 9 is approximated 
as constant, and v detection is considered possible for all zenith angles > 90°. 
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FIG. 22: Neutrino energy response for the up— going muon signal. The area under the curves is proportional to the signal of v induced 
muons (with threshold i? M = 1 GeV) assuming a power law neutrino flux with slope a. The thick (thin) curves are calculated for a slope 
a = 2 (a = 2.4); the solid (dashed) curves are calculated for a nadir angle 8 = 90° (9 = 0°) and negligible (maximum) v absorption in 
the Earth. The relative normalization of the curves is correct if the spectra have the same integrated flux above 1 TeV. The v induced 
signal is roughly proportional to the integrated neutrino flux in the approximate energy range [10 12 , 10 16 ] eV in case of a hard spectrum, 
or [10 , 10 ] eV for a soft spectrum. 
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FIG. 23: Neutrino to muon conversion efficiency. The different curves correspond to different choices for the muon detection threshold 
and different celestial declination of a source for a neutrino detector at the South Pole. 



